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ABSTRACT 

We present the results of 3D high resolution calculations of the last inspiral stages 
and the final coalescence of neutron star binary systems. The equations of hydrody- 
namics are solved using the smoothed particle hydrodynamics method with up to 10^ 
particles. Using Newtonian gravity, but adding the forces emerging from the emission 
of gravitational waves, we focus on the impact of microphysics on the dynamical evo- 
lution of the merger. Namely we use a new equation of state based on the relativistic 
mean field approach of Shen et al. ( 1998a, b). Neutrino emission of all flavours, the 
resulting cooling and the change in the electron fraction are accounted for with a de- 
tailed leakage scheme. 

The new equation of state is substantially stiffcr than the Lattimcr-Swesty equation of 
state that has been used in previous investigations. This leads the system to become 
dynamically unstable at a separation as large as 3.3 stellar radii, where the secu- 
lar orbital decay undergoes a transition towards a dynamical "plunge" of the binary 
components towards each other. As soon as the stars come into contact a Kelvin- 
Hclmholtz-instability forms at the interface of both stars. The peak temperatures are 
found in the vortex rolls that form during this process. We generally find slightly lower 
temperatures than previously found using the Lattimer-Swesty equation of state. 
The central object is surrounded by a very thick disk that shows cool equatorial flows: 
inflows in the inner parts of the disk and outflow further out. The cool inflows be- 
come shock heated in the innermost parts of the disk and lead to an outflow of hot 
material in the vertical direction. The temperatures in the disk have typical values 
of 3-4 MeV, lower than the temperatures found in previous investigations using the 
Lattimer-Swesty equation of state. These conditions allow for the existence of heavy 
nuclei even in the inner parts of the disk, we find typical mass fractions of ~ 0.1, 
which is enough for scattering off heavy nuclei to be the dominant source of neutrino 
opacity. 

The central object formed during the coalescence shows a rapid, differential rotation 
with periods of 2 ms. Although a final conclusion on this point is not possible 
from our basically Newtonian approach, we argue that the central object will remain 
stable without collapsing to a black hole, at least on the simulation time scale of 20 
ms, mainly stabilized by differential rotation. The massive, differentially rotating cen- 
tral object is expected to wind up initial magnetic fields to enormous field strengths of 
~ 10" G (Duncan & Thompson 1992) and may therefore have important implications 
for this event as a central engine of Gamma Ray Bursts. 

Key words: dense matter; hydrodynamics; neutrinos; gamma rays: bursts; stars: 
neutron; methods: numerical 



1 INTRODUCTION 

The final stage of the inspiral and the subsequent coales- 
cence of a neutron star binary system is among the prime 



candidates for ground based gravitational wave detection via 
the laser interferometers currently under construction such 
as LIGO (Abramovici et al., 1992), TAMA (Kuroda et al., 
1997), VIRGO (Bradaschia et al., 1990) and GEO (Danz- 
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mann, 1997). The neutron star merger scenario could also 
provide the energy reservoirs to power cosmological gamma 
ray bursts (GRB; Paczynski (1986); Eichler et al. (1989)) 
and is in addition one of the two most discussed production 
sites of the heavy, rapid neutron capture elements (Lattimcr 
& Schramm, 1974, 1976; Symbalisty & Schramm, 1982; Eich- 
ler et al., 1989; Rosswog et al., 1999; Preiburghaus et al., 
1999). Whether this scenario can fulfil all these promises or 
not, neutron star binary systems are known to exist and to 
inevitably coalesce (Taylor, 1994). This fact alone - apart 
from all further promises - justifies a careful study of the 
coalescence process and its possible implications. 
The coalescence is an intrinsically three-dimensional phe- 
nomenon and therefore analytical guidance is rare although 
very welcome and one has to resort to large scale compu- 
tations. Additional complications arise from the fact that 
there is almost no field of astrophysics that does not enter 
at some stage during the coalescence process: the last stages 
and the merger are certainly dominated by strong-field gen- 
eral relativistic gravity, the neutron star material follows the 
laws of hydrodynamics, particle physics enters via possible 
condensates of "exotic" matter in the high-density interiors 
of the neutron stars and the copiously produced neutrinos in 
the hot and dense neutron star debris, questions concerning 
element formation require detailed information on nuclear 
structure and reactions (often far from stability) to be in- 
cluded and also magnetic fields might play a decisive role 
since they may, via transport of angular momentum, de- 
termine whether and/or when the central, coalesced object 
collapses into a black hole. 

Due to this complexity current investigations follow one of 
two "orthogonal" lines: either ignoring microphysics, resort- 
ing to the simplest equations of state (EOS), polytropes, and 
thereby focusing on solving the complicated set of general 
relativistic fluid dynamics equations (or some approximation 
to it) or, along the other line, using Newtonian self-gravity 
of the fluid and investigating the influences of detailed mi- 
crophysics and relating the merger event to astrophysical 
phenomena. 

Many groups have performed 3D calculations of the merger 

scenario. The first work, using a polytropic EOS and Newto- 
nian gravity, was performed by Oohara & Nakamura (1997, 
and references therein) using Eulerian finite difference codes. 
Rasio & Shapiro (1992, 1994, 1995), Zhuge et al. (1994, 1996) 
and Davies et al. (1994) have used the smoothed particle 
hydrodynamics method (SPH) to further explore various as- 
pects of the scenario. 

The strong-field gravity aspect of the event was explored 
again by Oohara & Nakamura (1997, and references therein) 
using Post-Newtonian hydrodynamics in Eulerian formula- 
tion. Recently, Post-Newtonian calculations have also been 
performed using SPH by Ayal et al. (2001), Faber & Rasio 
(2000) and Faber et al. (2001). Several groups have worked 
on general relativistic formulations where the metric has 
been treated in the conformal flatness approximation (Wil- 
son & Mathews, 1995; Oohara & Nakamura, 1997; Baum- 
garte et al., 1997; Oechslin et al., 2001). Recently, much 
progress has been made towards a stable implementation of 
the fully relativistic equations of hydrodynamics (Shibata, 
1999; Shibata & Uryu, 2000). 

The line relying on (basically) Newtonian gravity and inves- 
tigating the microphysics of the event has been explored by 



Ruffert et al. (1996, 1997) and Ruffert & Janka (2001) using 
a grid code based on the piecewisc parabolic method (PPM) 
and by Rosswog et al. (1999, 2000) using SPH. 
In this paper we want to continue to explore the impact 
of the microphysics on the outcome of the coalescence. We 
use a new equation of state that overcomes previous restric- 
tions in the temperature and density. Since the high density 
regime of the nuclear physics is only poorly known to date 
this work will unfortunately not provide the final answer, 
but further broaden our understanding of the large effects 
that can be expected from the nuclear physics input. The 
calculations to be described below were performed on paral- 
lel computers allowing for an unprecedented hydrodynamic 
resolution. 

The paper is organized as follows. In Section 2 we describe 
the basic ingredients like the hydrodynamics or the equa- 
tion of state of our model. Section 3 describes the initial 

conditions of our calculations, the overall hydrodynamic evo- 
lution, and then focuses on aspects concerning the central 
object and the debris, referred to as "disk". The summary 
and a discussion of the results are provided in Section 4. 



2 MODEL 

2.1 Hydrodynamics 

To follow the dynamical evolution of the neutron star fluid 
we use a Lagrangian particle scheme, the so-called smoothed 
particle hydrodynamics method (SPH; Benz (1990); Mon- 
aghan (1992)). Since it is independent of any prescribed ge- 
ometry (e.g. grid) it is perfectly suited to handle the intrin- 
sically three-dimensional merger process. In addition, the 
Lagrangian nature of the scheme makes it easy to carefully 
track the evolution of interesting portions of the fluid (e.g. 
possible ejocta). Voids are treated in a natural way (i.e. no 
particles) and do not present any difficulty for the method, 
the interesting parts of the simulation do not have to be 
embedded in an artificial background medium like in other 
methods (e.g. PPM). The use of a background medium leads 
to further difficulties like emerging (artificial) shock waves 
at the stellar surfaces that have to be treated by additional 
remedies (Ruffert & Janka, 2001). If the corresponding pa- 
rameters are not chosen carefully, the artificial medium may 
also lead to a damping of (physical) oscillations (Font et al., 
2000). 

Since SPH is well-known we will only describe the ba- 
sic ingredients of our code. The basic set of equations may 
be found in Benz (1990). Rather than the "standard" form 
of artificial viscosity (AV; Monaghan & Gingold (1983)) we 
use a hybrid method that profits from two improvements: 
the viscosity parameters are time dependent (Morris & Mon- 
aghan, 1997) and have non-negligible values only in the pres- 
ence of shocks. Additionally, a switch is applied that sup- 
presses spurious forces in pure shear flows (Balsara, 1995). 
This hybrid method (Rosswog et al., 2000) has been shown 
to resolve shocks with the same accuracy as the standard for- 
mulation of AV, but exhibits much better behaviour in shear 
flows. In test calculations of differentially rotating stars the 
viscous time scales Tvisc,i = Vi /vi,visc obtained with the new 
scheme were two orders of magnitude longer than those from 
the standard form of AV. For details and test calculations 
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we refer to Rosswog et al. (2000). 

Since it is a widespread belief that SPH is always extremely 
viscous by nature, we want to measure the effective "alpha- 
viscosity", Oss, defined by Shakura & Sunyaev (1973) 



sCsH, 



(1) 



where u is the kinematic viscosity, Cs the sound velocity 

and H the disk thickness. In 3D the SPH bulk viscosity 
parameter, a (see e.g. Monaghan (1992)), can be related to 
the kinematic viscosity via (Murray (1995)) 



10 



Csh, 



(2) 



where h is the smoothing length. Due to the modifications 
in our hybrid AV-scheme we use 5i = Oi • fi in the above 
formula rather than the bulk viscosity parameter of particle 
i, Ui. Here fi is the Balsara- factor implemented in our code. 



fi 



IV ■ vl 



|V • v\i + |V X v\i + rics,i/hi 



(3) 



(4) 



and r] = 10~*. This factor suppresses AV in pure shear flows 

{fi ~> for |Vxt;| 3> |V-tT|) where it is unwanted, and fi^l 
in the case of pure shocks (|V x t7| <C | V • v|). Equating (1) 
and (2) one finds 

_ cti hi 

i.e. apart from the numerical value Oi it is the resolution of 
the disk that determines the effective alpha- viscosity ass ■ It 
should, however, be noted that this relation is only approxi- 
mate since the effects of AV and the a- viscosity prescription 

arc not necessarily exactly the same and due to the way the 
Balsara-factor is taken into account. To obtain numerical 
values we insert the average thickness of the central object 
and disk of run B (see below; Table 1) into equation (4). By 
averaging according to 



{a)ss = ^ 
where mi are the particle masses, we find 



(5) 
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For details of concerning "debris" and "central object" see 
below. These values are very low in comparison with the val- 
ues that are derived from observations of cataclysmic vari- 
ables, a ~ 0.1, and it is certainly possible that the physical 
viscosity of a hot neutron star debris disk is much higher 
than the values of the above determined numerical viscos- 
ity. 

The (Newtonian) forces of the self-gravity of the fluid arc 
efficiently calculated using a binary tree (Benz et al., 1990). 
Since the computationally most expensive part of the code 
is the evaluation of gravitational forces we implemented an 
integrator for the set of differential equations that only needs 
one force evaluation per time step. We decided to implement 
an Adams-Bashforth method which is third order accurate 
in time. This allows for a very accurate conservation of en- 
ergy and angular momentum even with relative large step 
sizes. In a calculation with ~ 100 000 particles these quan- 
tities are conserved to a few parts in 10* over ~ 10000 time 
steps. For reasons of comparison: other schemes (Ruffert & 



Janka, 2001) lose up to 10 % of the total angular momentum 
by rmrncrical artifacts despite excellent energy conservation. 
The amount of noise present in SPH depends, of course, on 
the particle distribution. For the worst stochastic 
particle distribution, it scales like Nn, where A'"jv is the 
number of neighbours a particle interacts with. Therefore 
while increasing the total particle number, N, one should 
also increase the number of neighbours in order to reduce 
the noise level. On the other hand one is not interested in 
too large an increase of the neighbour number and therefore 
the smoothing length, h, since it would compromise the spa- 
tial resolution and the numerical error; remember that the 
SPH-cquations are accurate to order O(h^). Therefore, A'"jv 
should increase slower than N. In the calculations presented 
here typical neighbour numbers are 100-120. 
Our whole code is parallelized for use on shared memory 
machines using OpenMP. For particle numbers above ~ 400 
000 we find an almost linear speedup for up to 120 proces- 
sors. 



2.2 Gravitational Wave Back-reaction 

The forces emerging from the emission of gravitational waves 
that drive the binary towards coalescence are treated in the 
point mass limit of the quadrupole approximation and are 
applied until the stars come into contact. Since each particle 
in a star is subject to the same back-reaction force the con- 
servation of the fluid circulation is guaranteed. For a further 
discussion of this approach we refer to Rosswog et al. (1999). 



2.3 Equation of State 

The equation of state (EOS) is one of the most crucial ingre- 
dients of a neutron star simulation and the microscopic be- 
haviour of matter decisively determines the overall, macro- 
scopic evolution of the system (compare, for example, the 
calculations using a polytrope versus those with the nuclear 
EOS of Lattimer & Swesty (1991) [LS-EOS] in Rosswog 
et al. (1999)). Since the stiffness of the EOS varies over a 
wide range as a function of density, temperature and com- 
position (see below) and the release of nuclear binding en- 
ergy when nucleons form nuclei has important dynamical 
consequences, a simple polytropic EOS is only a poor ap- 
proximation of the involved microphysics. All neutron star 
merger calculations that used a nuclear equation of state 
(Ruffert et al., 1996, 1997; Ruffert & Janka, 2001; Rosswog 
et al., 1999, 2000) relied on the Lattimer-Swesty-EOS. This 
EOS, however, has been designed for the use in supernova- 
calculations and therefore suffers from some deficiencies in 
our context: the electron fraction. Ye, is restricted to values 
above ~ 0.04, temperatures to values above ~ 10® K and 
densities above lO'^ gem ~^ . The electron fraction in a neu- 
tron star in /3-equilibrium, however, dips down to values of 
~ 0.01 near the neutron star surface (see below) and tem- 
peratures in old neutron stars are expected to be negligible 
with respect to typical nuclear energies; viscous heating dur- 
ing the inspiral is expected to increase the temperatures to 
only - 10* K (Lai, 1994). Therefore, once outside the EOS- 
range, these investigations either used extrapolation (Ruffert 
et al., 1996) or kept the EOS-quantities at the boundary val- 
ues (Rosswog et al., 1999). 
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In this work we use an EOS based on the tables of Shen 
et al. (1998a,b) that overcomes all above mentioned restric- 
tions, in the electron fraction as well as in temperature and 
density. Shen et al. follow a relativistic mean field (RMF) ap- 
proach which reproduces the basic features of more compli- 
cated relativistic Dirac-Briickner-Hartree-Fock calculations. 
The starting point is a relativistically covariant Lagrangian 
that contains, apart from the nucleon fields, the scalar o, 
and the uj and p vector mesons. For the a and uj meson fields 
non-linear terms are included, which are essential to repro- 
duce quantitative properties of nuclei (e.g. the compression 
modulus) correctly. The parameter set (TMl, Sugahara & 
Toki (1994)) is chosen in a way that experimental data of 
finite nuclei in the ground as well as in excited states are 
satisfactorily reproduced. With this parameter set the com- 
pression modulus, Kq, has a value of 281 MeV, whereas in 
the Lattimer-Swesty EOS in Ruffert et al. (1996), Rosswog 
et al. (1999) and all comparisons in this paper Kq = 180 
MeV has been used. Since many-body interactions at high 
densities are only poorly constrained to date, no effort is 
made to include more "exotic" physics such as hyperons, 
various mesons or quark matter in the high density regime. 
At densities above ~ 1/3 pnuc protons and neutrons form 
a homogeneous "nucleon fiuid", below this density matter 
may become inhomogeneous, i.e. the presence of nuclei may 
become energetically favorable. This phase is modelled us- 
ing the Thomas-Fermi approximation. Matter is assumed to 
consist of a mixture of nucleons, alpha-particles and spher- 
ical nuclei arranged on a lattice. The heavy nucleus (rep- 
resentative of a distribution of heavy nuclei) is assumed to 
be centered in a charge neutral cell consisting of vapor of 
the neutrons, protons and alpha particles. Above p ~ 10^^° 
g cm ~^ the nucleons are treated by the RMF theory, below 
this density they are assumed to form a Maxwell-Boltzmann 
gas. The alpha-particles are treated as a non-interacting 
Boltzmann gas, their occupied volume is accounted for in 
the calculation of the free energies. The density distribution 
in the Wigner-Seitz cell is parametrized and the free param- 
eters are determined by minimizing the free energy density 
with respect to the densities of all ingredients. 
Shen et al. only provide the baryonic part of the EOS. Their 
tables range from to 100 MeV in temperature, to ~ 0.56 
in Ye and 5.1 to 15.4 in \og{p) (p in cgs-units). We therefore 
add the contributions from photons and electron-positron 
pairs to the baryonic components. For the electron-positron 
pairs we use the code of Timmes & Arnett (1999). Apart 
from disregarding interactions, which is perfectly justified 
at the high densities of interest, the electron-positron-pairs 
are treated in full generality without any approximation. In 
the low density regime, p < 10^ g cm we extend the EOS 
with a gas consisting of neutrons, alpha-particles, electron- 
positron pairs and photons. The smooth transition is demon- 
strated in Figure 1. The EOS described above is tabulated 
with 32 entries in temperature, 72 in electron fraction and 
151 in density. We store the following quantities: total pres- 
sure, internal energy, difference in the nucleonic chemical 
potentials fi, sound velocity, entropy, mass fractions of pro- 
tons (xp), alpha particles (xa) and the heavy nucleus {xh), 
the proton to nucleon ratio, Z/A, and the nucleon number A 
of the heavy nucleus. The neutron mass fraction is obtained 

as Xn = 1 {Xp ~\~ Xa ~\~ Xh^ . 

In Fig. 2 we show a comparison between the pressures pro- 
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Figure 1. Pressure as a function of temperature for our EOS 
(baryonic part from Shen et al. -|- photons + electron-positron 
pairs, extension with a gas consisting of neutrons, alpha particles, 
electron-positron pairs and photons). The transition between the 
original EOS and the extension lies at p 10^ g/cm^. 
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Figure 2. Comparison of the pressures provided by the Lattimer- 
Swesty and the Shen et al. EOS as a function of density (cgs- 
units). 



vided by both the LS-EOS and the EOS of Shen et al. 
Below 10^^ gem ~^ the results are practically identical, at 
higher densities the new EOS is substantially stiffer than 
the LS-EOS. For rather low temperatures (T = 1 MeV) we 
find large differences in the density range from 10^^ to 10" 
gem which is the density regime that is crucial for the 
torus that forms around the central object, see below. For 
even higher densities the pressure curves for different tem- 
peratures converge indicating the diminishing relevance of 
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thermal pressure contributions. Again, the Shen et aL EOS 
yields considerably higher pressure values. 
The Shen-EOS also yields larger values for fi. Since the 
neutrino-less /3-equilibrium value of the electron chemical 
potential is given hy fis ~ fj- ^ Q, Q being the neutron pro- 
ton mass difference, and Ye is a monotonic function of p,e 
the Shen-EOS yields higher values for Ye than the LS-EOS. 



2.4 Neutrino Physics 

The emission of neutrinos provides a very efficient cooling 
mechanism for hot neutron star matter and the related elec- 
tron and positron captures change the electron fraction Ye 
which in turn determines the matter composition. We in- 
clude the most important neutrino reactions where we en- 
sure, via effective rates, that only the amount of neutrinos is 
produced that is actually able to leave the dense surround- 
ing material. Our scheme takes careful account of the energy 
dependence of the neutrino opacities by integrating over the 
neutrino distributions. 

For the emission processes we include electron captures (EC) 



e~ + p ^ n + Ve 

and positron captures (PC) 

e'^ + p + i>e 



(6) 



(7) 



which produce electron flavour neutrinos and the pair pro- 
ducing reactions, pair annihilation 



e +e'^ ^ Vi + Vi 
and plasmon decay 

7 — » + I'i- 



(8) 



(9) 



Ifcrc i'i/i'i refer to anti-/ncutriiio.s of all typos. 
For the opacities we include the dominant processes, neu- 
trino nucleon scattering 



Vi + {n,p} + {n,p}, 

coherent neutrino nucleus scattering 

and neutrino absorption by nucleons 

Ve + n ^ p -\- e~ 

I'c + p ^ n + e'^ . 



(10) 

(11) 

(12) 
(13) 



Note that this is the first time that the effects of scattering 
off heavy nuclei have been accounted for in a neutron star 
merger calculation. These arc important whenever nuclei arc 
present since the corresponding cross sections scale propor- 
tional to A^ and A « 80 (see below). For details of this 
transport scheme we refer to Rosswog et al. (2002), where 
the neutrino emission results are discussed in detail. 



3 RESULTS 

3.1 Initial Conditions 

3.1.1 Neutron Star Masses 

We focus in this investigation on equal mass systems with 1.4 
Mq per star since all well-determined neutron star masses 
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Figure 3. Profiles of matter density and electron fraction of the 
initial neutron star models (masses in solar units). 



from radio binary pulsars are distributed according to a nar- 
row distribution around 1.4 M© (Thorsett & Chakrabarti, 
1999). For neutron stars in X-ray binary systems substan- 
tially higher masses have been claimed: e.g. 1.8 ±0.2 M© for 
CygX-2 (Orosz & Kuulkers, 1999), ~ 1.9 M© for VelaX-1 
(J.H. van Kerkwijk & Zuiderwijk, 1995) and 1.8 ± 0.4 Mq 
for 4U 1700-37 (Heap & Corcoran, 1992). In order to get an 
upper limit we additionally explore the case of twice 2.0 M© 
(see Table 1). 

It should be noted that the system dynamics is extremely 
sensitive to oven small deviations from a mass ratio of unity. 
This sensitivity to the mass difference has been explored pre- 
viously, both using a polytropic (Rasio & Shapiro, 1994) and 
a realistic EOS (Rosswog et al., 2000; Ruffert & Janka, 2001) 
and will not be the subject of this investigation. 

3.1.2 Neutron Star Models 

We solve the one dimensional hydrostatic structure equa- 
tions together with the neutrino-less /3-equilibrium condi- 
tion to find the properties of the initial neutron stars. Ex- 
amples of such initial profiles are shown in Fig. 3. Note 
the enormous density decrease at the neutron star surface 
which is caused by the stiffness of the Shen-EOS. This effect 
is more pronounced for the higher mass neutron stars, low 
mass stars possess a thicker crust, as is visible from both 
density and Ve-profile. To illustrate this stiffness we show 
in Fig. 4 the quantity F = dln(p)/dln(p) obtained by finite 
differencing along a neutron star profile of 1.4 Mq. F rises 
from values close to 3 in the center to ~ 3.2 to drop sharply 
at the phase transition towards inhomogeneous matter and 
remain around 1.3 in the neutron star crust. 
To compare once again with the LS-EOS, we plot in Fig. 5 
the adiabatic exponent Fi, given by (Cox & Giuli, 1968) 



Fi = 



du\- 
dT) 



(14) 



where T is the temperature and u the specific internal en- 
ergy, for both the Shen-EOS and the LS-EOS in the relevant 
low- Ye regime. Ye = 0.1. In the case of the Shen-EOS the 
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Table 1. Summary of the different runs, ao : initial separation; f : neutrino physics; Tgj^ : simulated duration; M1/M2 : masses in solar 
units; # part.: total particle number 



run spin Mi M2 # part. ao [km] v Tgim [ms] 



A 


corot. 


1.4 


1.4 


207,918 


48 


no 


10.7 


B 


corot. 


1.4 


1.4 


1,005,582 


48 


no 


10.8 


C 


irrot. 


1.4 


1.4 


383,470 


48 


yes 


18.3 


D 


corot. 


1.4 


1.4 


207,918 


48 


yes 


20.2 


E 


irrot. 


2.0 


2.0 


750,000 


48 


yes 


12.2 



_ 1 1 1 1 1 1 1 1 1 1 1 1 1 
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derivatives in (14) have to be obtained by finite differenc- 
ing, in the LS-EOS case they are available from analytical 
expressions. For the cold case, T = 0.5 MeV, Fi rises from 
values of ~ 1.5 at log(p) = 10 to values of ~ 3.3 for the 
Shen-EOS and ~ 3.1 for the LS-EOS. Note the sharp peak 
at log(p) m 14 which is present in both EOSs. In the high 
temperature case, T = 10 MeV, the maximum values axe 
slightly lower and the peaks are smeared out. This variation 
of Fi with the physical conditions demonstrates again the 
difficulty of approximating a nuclear EOS by a polytrope of 
fixed adiabatic exponent. 

The electron fraction, Ye, decreases from central values of 
around 0.1 (higher values arc encountered for more massive 
stars) to values as low as ~ 0.01 and increases steeply again 
towards the surface (Fig. 3). Note that only a tiny amount 
of material, ~ 1% of the star's mass, is located in this region 
of increasing Ye. 



Figure 4. F = dln(p)/dln(p) along a neutron star profile of 1.4 
M0. 



0-© Shen et al., Ye= 0.1, T= 0.5 MeV 
— Shenetal.,Ye=0.1,T=10MeV 
1 LS-EOS, Ye= 0.1, T= 0.5 MeV 
LS-EOS, Ye= 0.1, T= 10 MeV 
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Figure 5. Shown is the variation of Fi with density (cgs-units) 
and temperature for a fixed electron fraction. Ye = 0.1, for 
both our equation of state (Shen et al. plus extension) and the 
Lattimer-S westy-EOS . 



3.1.3 Initial Spins and Separation 

Since the time during which the neutron stars can tidally in- 
teract is extremely short it would need an implausibly high 
neutron star viscosity to lead to a tidal locking of the binary 
components during the inspiral phase (Bildsten & Cutler, 
1992; Kochanek, 1992). The spins at the moment of merger 
will be negligible with respect to the orbital angular mo- 
mentum and therefore the most realistic spin configuration 
is close to the irrotational case. Thus, the most interesting 
configuration to be explored is the case without initial neu- 
tron star spins. We also investigate corotating systems for 
reasons of completeness and since it is straightforward to 
construct equilibrium configurations by damping out veloc- 
ities in the corotating frame. For details of this procedure 
we refer to Rosswog et al. (1999). The case of neutron stars 
spinning in the direction opposite to the orbit has been ex- 
plored in previous work (RufFert et al., 1996; Rosswog et al., 
1999) and shall not be further discussed here. 
The neutron stars with the corresponding spins arc then set 
to Keplerian orbits and are provided with radial velocities 
corresponding to 



ao 



64 M^fx 

5 a?. ' 



(15) 



where ao is the initial separation, M the total and the 
reduced mass (e.g. Shapiro & Teukolsky (1983)). For fur- 
ther details we refer to Rosswog et al. (1999). Starting with 
spherical stars is obviously not an equilibrium configuration 
and will therefore result in oscillations of the stars. We have 
tried to reduce these oscillations by stretching the stars ac- 
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Figure 6. Mass distribution in the orbital plane. The labels indicate contours of log(p), p is in gcm~^. The left column corresponds to 
a corotating system (more than 10^ particles, run B) the right one to a system without initial neutron star spins (more than 700 000 
particles, run E). The orbital motion is counter-clockwise. 



Figure 7. Left column: velocity field and density contours in the orbital plane in the corotating frame for an initially corotating 
configuration (> 1 000 000 particles, run B). Shown is the central object only (down to log(p= 13.5, p in g cm~*). Right column: 
corresponding temperatures. 



cording to the ellipsoidal approximation of Lai et al. (1994a) 
for polytropcs of F = 3 (compare Fig. 4) and correcting for 
the finite size effects in the corresponding orbital frequency. 
However, this did not improve the calculations, we therefore 
started the calculation with spherical stars and Keplerian 
orbital frequency. At the chosen initial separation (see be- 
low) the tidal deformations of the stars are very small, the 
height of the tidal bulge being h w (ii*/oo)^/J* w 0.03 R», 
where R» is the neutron star radius. 

The initial separation has to be determined as a tradeoff 
between computational resources and physically reasonable 
initial conditions. Binary systems can become dynamically 
unstable (i.e. they coalesce on a dynamical time scale) due to 
entirely Newtonian tidal effects (Chandrasekhar, 1975; Tas- 
soul, 1975). This means that orbital decay changes abruptly 
from the secular orbital decay driven by the emission of grav- 
itational waves to a rapid plunging of both components to- 
wards each other on just the orbital time scale. The reason 
for this instability is the steeping of the effective interac- 
tion potential between both components due to tidal effects. 
Tidal effects increase with the incompressibility of the stel- 
lar fluid (since the stars are less centrally condensed) and 
therefore the onset of this dynamical instability is a very 
sensitive function of the stiffness of the EOS, setting in at 
larger separations for a stiffer EOS. This instability has been 
studied extensively in Newtonian gravity using an ellipsoidal 
approximation to the neutron star shape (Lai et al. (1994a) 
and references therein). Rclativistic effects shift the inner- 
most stable circular orbit, -Risco, to larger binary separa- 
tions (Baumgarte et al., 1998a,b). 

We determined the binary separation where the system be- 
comes dynamically unstable experimentally. To this end a 
binary configuration was relaxed in the mutual tidal field, 
set to a corotating orbit and evolved for several orbits with- 
out radial velocity and back-reaction force. This experiment 
again underlined the enormous stiffness of the EOS which 
translated in a (for a Newtonian calculation) very large sep- 
aration for the ISCO. We found a corotating binary to be 
stable for an initial separation of Rdyn = 49.5 km, corre- 
sponding to JZisco/^* 3.3. We therefore chose the initial 
separation for the simulation start as ao = 48 km, which we 
regard to be an acceptable compromise between computa- 
tional effort and physical appropriateness. 
The simulations performed are summarized in Table 1. We 
use particle numbers up to 10® which translates into smooth- 
ing lengths of w 0.38 km in the initial neutron stars. 

3.2 Overall Hydrodynamic Evolution 

Since we start our simulation just inside the last stable 
orbit the neutron stars approach each other very quickly 
and merge within roughly one orbital revolution. Being con- 



stantly kept back by a slowly receding centrifugal barrier 
the merger itself proceeds very subsonically (keep in mind 
that typical sound velocities are ~ 0.4 c). 
The corotating systems (runs A, B, D) merge within approx- 
imately one orbital period. Prior to merger only a tiny lag 
angle develops between the axes of the neutron stars (first 
panel, left column Fig. 6). Note that such a lag angle devel- 
ops even in absence of viscosity since the system is not able 
to adapt fast enough to the rapidly changing tidal potential 
(Lai et al., 1994b). Immediately after contact mass shed- 
ding via the outer Lagrangian points sets in which results 
in thick, puffed up spiral arms (Fig. 6), that subsequently 
wrap around the central object to form a disk. The spiral 
arms show an appreciable lateral expansion and soon engulf 
the central object and the innermost high-density parts of 
the disk (last panel, left column in Fig.6). 
The irrotational configurations merge slightly quicker, the 
system with twice 1.4 M© (run C) after around three- 
quarters of an orbit and the system with 2.0 Mq (run E; 
Fig.6) after half an orbit, due to the lower total angular 
momentum. This fgister inspiral leads to noticeably larger 
lag angles before contact. The systems start immediately to 
shed hot material with rather low density from the interac- 
tion region. This is closely followed by mass shedding via 
the outer Lagrangian points. The matter that is shed by 
the latter mechanism whips through the previously ejected 
material forming a spiral shock which heats up the mate- 
rial. The spiral pattern is much less pronounced than in the 
corotating case and becomes washed out within a few mil- 
liseconds by lateral expansion and a rapidly expanding disk 
is left behind. The material within this disk follows eccen- 
tric trajectories and once the outward motion is reversed and 
matter starts falling back towards the hot compact remnant 
in the centre, the inner parts of this torus-like structure be- 
come compressed and heat up. Note the enormous expansion 
of the debris material, in all cases the typical diameter of the 
mass distribution increases from ~ 100 km to ~ 1000 km in 
just a few milliseconds. 

One may speculate that the stronger general relativistic 
gravity would lead to a more compact post-merger configu- 
ration with a more massive central object and a less massive 
torus around it. For polytropic calculations this tendency is 
indeed visible (Shibata & Uryu, 2001). 

In Table 2 we show the (baryonic) masses contained in the 
central object and the debris material. Both regimes are 
discriminated on grounds of a log(p)-r-plot. In several cases 
the transition between central object and debris is rather 
smooth so that it is difficult to determine the exact mass 
cut. The resulting uncertainties are of the order 0.02 Mq. 
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Figure 8. Left column: velocity field and density contours in the orbital plane in the corotating frame for a system without spins 
(~ 400 000 particles, run C). Shown is the central object only (down to log(p= 13.5, p in g cm~^). Right column: corresponding 
temperatures. 



Table 2. Final rest mass distribution, Mco refers to the mass of 
the central object and Mjeb the mass of the debris material, a 
denotes the stability parameter Jc/GM^^, where J is the angular 
momentum contained in the central object. 



run 


Mco 


Mdeb 


a 


A 


2.33 


0.47 


0.48 


B 


2.35 


0.45 


0.49 


C 


2.55 


0.25 


0.64 


D 


2.33 


0.47 


0.48 


F 


3.71 


0.55 


0.55 



0.06 
< to > 

0.04 



pos. X-axis, z= km 
pos. X-axis, z= 7.5 km 
diag. (pos. X, neg. y), z= km 
diag. (pos. X, neg. y), z= 7.5 km 



3.3 Central Object 

The central objects of all performed calculations do not 
become axisymmetric on the simulation time scale (Figs. 
6,7,8) and will therefore continue to emit gravitational radi- 
ation. This result is in qualitative agreement with those from 
uniformly rotating stiff polytropes (F > Fcrit ~ 2.3; Tas- 
soul (2001)). Note that contrary to previous investigations 
(Davies et al., 1994; Rasio & Shapiro, 1994) we never end up 
with rigidly rotating central objects. This is an effect from 
the lower viscosity due to our new artificial viscosity scheme 
and the high resolution of the present calculations; the ef- 
fective alpha- viscosities in the central object are ~ 10~®, see 
Section 2.1. 

It has been pointed out previously by various authors (Ruf- 
fert et al., 1996; Rosswog et al., 1999; Rasio & Shapiro, 1999; 
Faber et al., 2001) that a vortex sheet forms between the 
stars as soon as they come into contact. This is most pro- 
nounced for the irrotational case, since, in the frame corotat- 
ing with the binary, the stars seem to be spinning in opposite 
directions and therefore a discontinuous velocity field is en- 
countered when the interface between both stars is crossed. 
Such a vortex sheet is known to be Kelvin-Helmholtz un- 
stable on all wave lengths, with the shortest modes grow- 
ing fastest. As has been pointed out previously, the shortest 
growing mode is determined by the smallest, numerically re- 
solvable length scale. For a further discussion of this point 
we refer to the previously mentioned literature. 
The velocity field for the initially corotating system is shown 
in column one in Fig. 7. Along the vortex sheet two vortices 
form which remain well-separated and which are not dissi- 
pated until the end of the simulation. In Fig. 8 the velocity 
fields in a frame corotating with the binary system are shown 
(orbital plane, central object only) for the spinless system. 
The vortex sheet is clearly visible in the first panels of Figs. 
7 and 8. The Kelvin-Helmholtz instability creates a set of 
vortices along the interaction interface. These merge during 
the further evolution leaving behind a rapidly, differentially 
rotating merger remnant. As will be discussed below, this 
has important consequences for the stability of the central 
object. 

To further illustrate the diflferential rotation of the merged 



5 10 15 20 25 

fcyl M 

Figure 9. Differential rotation of the merged remnant I: snapshot 
of the angular frequencies (displayed in code units of 1.984 • 10® 

1/s), of the most realistic initial configuration, run C, at t= 7.686 
ms (compare Fig. 8) along the positive x-axis (for two different 
heights) and along a diagonal (positive x, negative y; again two 
heights) . The increase of w towards the origin corresponds to the 
central vortex (see Fig. 8). 

remnant we show in Fig. 9 the SPH-smoothed values of the 

angular frequency 

(^^)^r) = J2^^W{\r-r-\,hj), (16) 

J 

where r is the point of interest (see below), Vtan,j is the cylin- 
drical tangential velocity of particle j, Vcyij = -s/x^ + f | > 
nij and pj are mass and density and W is the used spline- 
based SPH-kcrncl (e.g. Monaghan (1992)). We only show 
the case of no initial spins (run C), since it is the most rel- 
evant one (the corotating case is apart from the persistent 
vortices closer to rigid rotation). The chosen time, t= 7.686 
ms, corresponds to the last panels, column one and two in 
Fig. 8. We show (ui) along the positive x-axis, both at the 
height of z= and z= 7.5 km, and along a diagonal given by 
positive x-values and negative y-values, since this line is an 
approximate symmetry-axis at this stage (again for heights 
of z— and z— 7.5 km), (a;) is strongly peaked at the centre, 
where the central vortex is located and decays with increas- 
ing distance from the origin depending on the direction and 
height above the orbital plane. The corresponding periods, 
(T) = 2tv/{lj) arc shown in Fig. 10. The central 10 km of all 
chosen lines have periods below 1 ms. 

The highest temperatures of the merged configuration 
are found in the above mentioned vortices, see right columns 
Figs. 7 and 8. The maximum temperatures inside the cen- 
tral objects arc shown in Fig. 11. Due to the high degener- 
acy of the neutron star matter oscillations due to imperfect 
initial conditions and numerical noise lead to temperatures 
of a few MeV, however only in the dense, central regions 
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Figure 11. Maximum temperatures (in MeV) in the merged con- 
figurations as a function of time (ms). 

where these finite temperatures have no physical efi'ect (no 
thermal pressure contribution at these densities; matter is 
opaque to neutrinos) . In the case of corotation the tempera- 
ture hardly rises above the level during inspiral. Due to the 
larger shear motion irrotational initial conditions yield sub- 
stantially higher maximum temperatures; up to ~ 30 MeV 
in the case of 1.4 Mq and ~ 45 MeV for our extreme case 
of 2.0 Mq stars. These temperatures are lower than those 
reported by Ruffert et al. (1996). This may be an effect of 
the different masses used in both investigations (1.4 versus 
1.6 Mq), the different EOS (compare also Ruffert & Janka 
(2001), Fig. 15) and different numerical viscosities in both 
codes. 

We find the masses of the central objects to be roughly 
0.1 Mq lower, sec Table 2, than those of previous calcula- 
tions using the softer Lattimer-Swesty EOS (Rosswog et al., 
1999). 

The question whether/when the central object collapses to 



a black hole is relevant for various reasons, e.g. to decide 
on possible mechanisms to produce Gamma Ray Bursts 
(GRBs). Most state of the art nuclear equations of state can 
support cold, non-rotating neutron stars in /3-equilibrium of 
~ 2.2 Mq gravitational mass (Akmal et al., 1998). This cor- 
responds approximately to the baryonic masses we find here 
for the central objects. Recent relativistic mean field equa- 
tions of state are able to stabilize even heavier neutron stars 
(between 2.45 and 3.26 Mq, (Chung et al., 2001)). There- 
fore even cold, non-rotating configurations of the masses 
found here may be supported by most of these EOSs. Addi- 
tional stabilizing effects come from the finite, thermal con- 
tributions to the pressure. If the matter further contains 
non-Ieptonic, negative charges such as hyperons, d or 
s quarks, then trapped neutrinos lead to an additional in- 
crease of the maximum mass (Prakash et al., 1995). The 
initial magnetic field of neutron stars, So ~ 10^^ G, is ex- 
pected to be amplified in the differentially rotating remnant 
to enormous field strengths, ~ lO^'^ G (Duncan & Thomp- 
son, 1992), which, again, leads to a substantial increase of 
the maximum mass (Cardall et al., 2001). Stark & Piran 
(1985) have analyzed the stability of a rotating polytrope 
against the collapse towards a black hole. For the investi- 
gated case of rigid rotation they found that rotation can 
stabilize the fluid against collapse if the stability parame- 
ter o = Jc/GM'^ is larger than w 1. The stability param- 
eters a of the central objects in our calculations are ~ 0.5 
and arc shown for completeness in Table 2. We expect the 
most important stabilizing effect to come from the differ- 
ential rotation of the merger remnant. The enormous in- 
crease in the maximum mass by differential rotation has 
boon demonstrated for the case of white dwarfs by Ostriker 
& Bodenheimer (1968) who constructed configurations up 
to 4.1 Mq. Recent general relativistic treatments of differ- 
entially rotating neutron stars (Baumgarte et al., 2000) find 
that even modest degrees of differential rotation allow for 
significantly higher maximum masses than non or uniformly 
rotating configurations. 

It may be somewhat bold to speculate from an essentially 
Newtonian calculation about the stability of the remnant, 
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but from the above line of argument we are tempted to con- 
clude that the central object remains stable at least on the 
simulation time scale, ~ 20 ms. Due to our ignorance of 
the high-density equation of state, it cannot be ruled out 
that the end product of the coalescence is a stable super- 
massive neutron star of ~ 2.8 Mq. It seems, however, more 
likely that the central object is only temporarily stabilized 
and once the stabilizing effects weaken (e.g. neutrinos have 
diffused out after ~ 10 s, magnetic braking has damped 
out differential rotation) collapse to a black hole will occur. 
The time scale until collapse is difficult to determine, since 
it depends sensitively on poorly known physics and on the 
specific system parameters. Wc expect the most important 
effect to come from the differential rotation Ostrikcr & Bo- 
denheimer (1968); Baumgarte et al. (2000) and therefore the 
collapse time scale to be set by the time it takes the remnant 
to reach uniform rotation. Assuming the dominant effect to 
come from magnetic dipole radiation (the viscous time scale 
is estimated to be ~ 10^ s (Shapiro, 2000)), the time till 
collapse is estimated as 

18c-'*M 

-.^2 ( M \ Ao^'^gN Vl5km\ V3000s-^\'',_, 

where M, B and R are mass, magnetic field and radius of 
the central object. 



3.4 Disk 

Immediately after contact the two merged stars try to get 
rid of excess angular momentum by shedding mass. Wc refer 
to this debris material as 'disk'. In all cases, a few millisec- 
onds after contact a thick disk, with typical heights roughly 
comparable to the disk radius, has formed around the cen- 
tral object (see Figs. 6, 14 and 15). The masses contained 
in the debris material are, depending on the initial spin, be- 
tween 0.25 for run C and 0.55 M© for our extreme case, run 
E. 

Initially, as the material is being shed by the central object, 
the debris material is cooled via the very fast expansion. 
Part of this material is later shock-heated as subsequently 
ejected material runs into it. This is most pronounced in the 
case of initial corotation. Here the spiral arms expand later- 
ally due to the released nuclear energy from recombination 
of nucleons into nuclei and due to the change in the adiabatic 
exponent and hit supersonically. This heating mechanism is 
illustrated for a run A in Fig. 13. 

Fig. 14 shows disk properties in the orbital plane for the 
most realistic, initially non-spinning, configuration. The 
shedded mass follows eccentric trajectories, resulting in a 
fast initial expansion and cooling phase during the outward 
motion (see panel one, left column in Fig. 14) and a sub- 
sequent compression when the motion is reversed and the 
inner parts start to fall back towards the central object (see 
panel one, column two in Fig. 14), against material that is 
still being shed and moving outwards. 

To illustrate the interplay between the fluid motion, the tem- 
perature and the vertical disk structure, we show in Fig. 15 
the velocity fields together with density contours (left col- 
umn) for run B (twice 1.4 Mq, corotation), run C (twice 1.4 



M0, no initial spins) and our extreme case, run E (twice 2.0 
M0, no initial spins). The corresponding temperatures are 
shown in the right column. The disk shows violent convec- 
tive motion, that can be subdivided into roughly four flow 
regimes: (i) the innermost, shock- heated 'torus', that results 
from the collision of the hot matter shed from the central 
object and the cool (T < 1 MeV) equatorial inflow; (ii) the 
cool, equatorial inflow phase; (iii) a rarefaction regime that 
separates the inflow from the outflow and (vi) the cooling 
outflow from the disk. As the inflow hits the torus it drives a 
hot, vertical outflow which leads to the thick (H~R), puffed 
up disk. Only the innermost, high density material of the 
debris has temperatures in excess of ~ 4 McV, the shocked 
torus is at around 3 MeV, the hot, high-altitude outflow 
with densities below 10^ gcm^"* has temperatures between 1 
and 2 MeV, the outer parts of the flow close to the equato- 
rial plane are substantially cooler (see also Fig. 14, column 
two, panel two). In Fig. 16 we show a blow-up of the veloc- 
ity flelds (run B and E) in the regime where the cool inflow 
hits the inner parts of the disk and produces the "butterfly- 
shaped" temperature distribution. The temperatures found 
are slightly lower than those reported by RufTert & Janka 
(2001) (up to 10 MeV). Wc suspect the stiffness of the Shcn- 
EOS in the range between 10^^ and 10^'* g cm~^, see Fig. 2, 
to be partly responsible for this fact. 

Despite the high temperatures we still find a substantial 
meiss fraction of heavy nuclei in the inner parts of the disk, 
which dominate the neutrino emission. Wc find roughly 10 
% of the disk mass in the form of heavy nuclei (third panels 
in columns one and two of Fig. 14). The Shen-EOS yields 
in the inner torus regime an average nucleus with yl ~ 80 
and Z ~ 0.3. We find for the ratio of the mean free path 
due to scattering off heavy rmclei, Aa, and off free rmcle- 
ons. An, Aa/A„ = n„(j„/{nAcyA) ~ 4(1 - Xh)/{XhA) « 0.5. 
Therefore the scattering mean free path in these regions is 
dominated by heavy nuclei. These have not been accounted 
for in previous investigations, which assumed the torus to 
be fully photo-disintegrated (Ruffert et al., 1996). 
In Fig. 17 we plot the logarithm of the angular velocities 
Eis a function of the cylindrical radius for runs B, C and E 
and compare it to the corresponding Keplerian angular ve- 
locities. For runs B and C the innermost parts of the disk, 
up to ^ 100 km, arc roughly Keplerian. Further out, up to 
~ 300 km the disk rotates sub-Keplerian and the outermost, 
outflowing parts rotate faster than the corresponding Kep- 
lerian value. The disk of run E is roughly Keplerian out to 
~ 200 km, drops below the Kepler value between ~ 200 and 
~ 450 km and is above it at larger radii. 



4 SUMMARY AND DISCUSSION 

We have performed high-resolution 3D calculations of the 

last moments of the inspiral and the final coalescence of a 
neutron star binary system. Our main motivation was to 
explore the impact of the, sometimes poorly known, micro- 
physics on the merger scenario. 

The equations of hydrodynamics were solved using the 
smoothed particle hydrodynamics method with particle 
numbers of more than 10® to resolve the details of the 
thermodynamic and nuclear evolution. The initial neutron 
stars are resolved with smoothing lengths of w 0.38 km. We 
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Figure 13. Disk heating mechanism: temperatures of the central region shortly after the merger (t= 5.225 ms, corotation, ~ 200000 
particles, run A). The disk is shock-heated when the spiral arms hit supersonically due to their lateral expansion (green). The two hot 
spots (dark blue) in the central object stem from the vortices which formed during the Kelvin-Helmholtz instability during the merger 
(see text). 



Figure 14. Physical properties of the disk (equatorial plane) that forms around the central object (no initial spins, run C). Shown are 
snapshots t= 7.434 ms (left column) and t= 11.214 ms (right column). 



have measured the effective alpha-viscosity of our code and 
find for the well-resolved calculations values of < a >ss~ 
6 • 10"''. The self-gravity of the fluid was treated in a New- 
tonian fashion, but forces from the emission of gravitational 
waves were added. 

In these calculations the temperature-dependent, micro- 
scopic behaviour of the matter is modelled using a new nu- 
clear equation of state whose baryonic part is described in 
the relativistic mean field approach (Shen et al., 1998a,b). 
We have added the contributions of electron-positron pairs 
(treated without any approximation) and photons to the 
baryonic component. In the low density regime the equation 
of state has been extended with a gas consisting of neu- 
trons, alpha particles, electron-positron pairs and photons. 
This new EOS covers the whole parameter space in density, 
temperature and electron fraction that is relevant to the 
neutron star merger scenario and therefore overcomes the 
restrictions inherent to previous calculations using a differ- 
ent equation of state (Ruffert et al., 1996; Rosswog et al., 
1999). 

We have also included the change in the electron fraction. 
Ye, and internal energy due to the emission of neutrinos by 
means of a detailed multi-fiavour neutrino treatment that 
accounts carefully for the energy dependence of the neutrino 
reactions. In contrast to previous calculations scattering off 
heavy nuclei was included as a source of opacity. 
The new equation of state is substantially stiffer than the 
previously used Lattimer-Swesty EOS, with adiabatic expo- 
nents reaching values well above 3 within a 1.4 Mq star. This 
has several important implications. First, since the stars are 
less centrally condensed the binary system becomes dynam- 
ically unstable at a larger separation. Since the inspiral be- 
haviour changes qualitatively at this point from a quasi- 
stationary, quasi-periodic motion to a dynamic "plunge" to- 



wards each other the gravitational wave signal will change 
correspondingly at this point. For a corotating system of 1.4 
M0 we find the instability to set in at a separation of ~ 
49 km and an orbital frequency of ~ 280 Hz corresponding 
to « 560 Hz in the gravitational wave signal, well within 
the range accessible to GEO600 and LIGO (Schutz & Ricci, 
2001). 

Once the neutron stars come into contact a Kelvin- 
Helmholtz unstable vortex sheet forms between both stars 
leading to the formation of vortex rolls. In the case of an ini- 
tially corotating system two such vortices form which remain 
present and well-separated until the end of the calculation. 
In the more important irrotational case, several vortices 
form which merge within an orbital time scale leaving be- 
hind a differentially rotating central object. It is within these 
vortex roles that the highest temperatures of the merged 
configuration are found, reaching peak values of 30 MeV for 
the case of an irrotational system, while in the corotating 
case only a moderate temperature increase is observed. In 
a test calculation using two stars of 2.0 Mq temperatures 
of up to 45 MeV are reached. We find generally comparable 
but slightly lower temperatures than previous investigations 
that used the Lattimer-Swesty EOS. 

The disks have a vertical height comparable to their radial 
extension. We find the vertical disk structure to be deter- 
mined by an equatorial inflow of cool material that becomes 
shock heated in the innermost, dense disk regime. This goes 
along with an outflow of hot material in the vertical direc- 
tion. In the dense inner torus around the central object we 
find temperatures of 3-4 MeV as compared to 5-10 MeV 
in the Lattimer-Sesty case. Despite the high temperatures 
we find a substantial amount (a mass fraction of ~ 10%) 
of heavy nuclei present even in the inner parts of the torus. 
This is enough to dominate the neutrino scattering opacities 
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Figure 15. Vertical disk structure; density contours and velocity fields of runs B, C and E at around 10.7 ms are shown in the left 
column; the right column shows the temperature distribution (densities above 10^^ gcm^ are cut out to enhance the contrasts). 




Figure 16. Blow-up of the velocity fields of the inner disk regions of run B and E (compare to Fig. 15). 




r cyl [km] 



Figure 17. Plotted is the logarithm of the angular velocities (in 
code units of 1.984 • 10^ 1/s) along the y-axis as a function of the 
cylindrical radius for runs B (circles; t= 10.685 ms), C (squares; 
t= 10.710 ms) and E (triangles; t= 10.710 ms). This is compared 
to the corresponding Keplerian angular velocities. 



which therefore play a crucial role for the neutrino transport. 
This has important consequences for the enormously tem- 
perature sensitive neutrino reactions and will be discussed 
in a future paper. 

The masses of the central objects arc roughly 0.1 M© lower 
than in previous calculations where we used the Lattimer- 
Swesty-EOS, leading to masses of the central object that 
are comparable to the maximum masses that most recent 



nuclear equations of state can support for cold, non-rotating 
configurations. Since several stabilizing efi^ects like thermal 
pressure contributions, rotation and, most likely, magnetic 
fields axe active during the coalescence, we suspect the cen- 
tral object will remain stable on at least the simulation 
time scale of 20 ms. Probably the most important stabil- 
isation comes from differential rotation. For the most prob- 
able case, where the neutron star spins are negligible with 
respect to the orbital motion, the maximum density of the 
central object barely reaches the initial central density of 
a single neutron star although its mass is roughly one so- 
lar mass higher. We estimate that the central objects could 
be stabilized against collapse for as long as ~ 10^ s. It has 
to be stressed, however, that none of the equations of state 
used so far contains particles more exotic than nucleons. 
The appearance of more exotic matter in the deep cores of 
neutron stars is still subject to large uncertainties due to 
poorly constrained many-body interactions at highest den- 
sities. Therefore final conclusions on this point cannot be 
drawn. 

If the central object indeed remains stable magnetic fields 
may wind up to reach ~ 10^^ G for a typical rotation pe- 
riod of 2 ms (Duncan & Thompson, 1992) and may therefore 
provide the conditions for magnetically-powered gamma ray 
bursts. 
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